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DIRECT COUPLING OF FLUIDS AND STRUCTURES 
FOR AEROSPACE APPLICATIONS 


I. INTRODUCTION 

In recent years, significant advances have been made for single disciplines in both 
computational fluid dynamics (CFD) using finite-difference approaches 1 and computa- 
tional structural dynamics (CSD) using finite-element methods (see Chapter I of Ref. 2). 
For aerospace vehicles, structures are dominated by internal discontinuous members such 
as spars, ribs, panels and bulkheads. The finite-element (FE) method, which is fundamen- 
tally based on discretization, has proven to be computationally efficient to solve aerospace 
structures problems. The external aerodynamics of aerospace vehicles is characterized 
by field discontinuities such as shock waves and flow separations. Finite-difference (FD) 
computational methods have proven to be efficient to solve such problems. 

However, only a limited amount of work has been completed in coupling these two 
disciplines for multidisciplinary applications. The prime reason is the lack of computational 
power for combining the two major computational fields. The development of a new 
generation of parallel computers can possibly alleviate the restriction of computational 
power. 

In addition, such coupling procedures result in an increased level of complication. 
Therefore, aeroelastic analysis has been mostly performed by coupling advanced computa- 
tional fluid dynamics (CFD) methods with simple structural modal equations or advanced 
computational structural dynamics (CSD) methods with simple flow solutions. However, 
these approaches can be less accurate for the aeroelastic analysis of practical problems 
such as a full aircraft configuration in the transonic regime. Moreover aeroelastic prob- 
lems of aerospace vehicles axe often dominated by large structural deformations and high 
flow nonlineaxities. It is necessary to develop a fully coupled procedure utilizing advanced 
computational methods for both disciplines. 

The objective of this research is to develop computationally efficient methods for solv- 
ing fluid-structural interaction problems by directly coupling finite difference Euler/Navier- 
Stokes equations for fluids and finite element dynamics equations for structures on parallel 
computers. This capability will significantly impact many aerospace projects of national 
importance such as Advanced Subsonic Civil Transport (ASCT), where the structural sta- 
bility margin becomes very critical at the transonic region. This research effort will have 
direct impact on the High Performance Computing and Communication (HPCC) Program 
of NASA in the area of parallel computing. 
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II. PREVIOUS STATUS 


A multidisciplinary code for computing unsteady flows and aeroelastic responses of 
aerospace vehicles, ENSAERO, has been developed on serial supercomputers at the Com- 
putational Aerosciences Branch of the NASA Ames Research Center. 3 This multidisci- 
plinary code computes unsteady aerodynamic responses of aircraft using the Euler/Navier- 
Stokes equations. An aeroelastic shape- conforming moving grid is used to include the effect 
of structural deformations on unsteady flows. This code is designed in a modular fashion 
to adopt several different numerical schemes suitable for accurate aeroelastic computa- 
tions. The basic coding of ENSAERO can accommodate zonal grid techniques for efficient 
modeling of full aircraft. 

An early version of ENSAERO 4 has been successfully applied in computing aeroelas- 
tic responses of a rectangular wing by using the Euler equations for fluids and the modal 
equations for structures. The result demonstrates that the code can accurately predict the 
flutter dynamic pressure of a rectangular wing. The code was extended to compute aeroe- 
lastic responses using the Navier-Stokes equations for fluids. 5 Later, it was updated by 
utilizing an upwind algorithm, and the code has been applied to fighter wings undergoing 
unsteady motions 6,7 at moderately large angles of attack. This code also has a capability 
of modeling moving control surfaces. 8 Furthermore, ENSAERO has demonstrated the ca- 
pability to simulate transonic flows on wing-body configurations using the Navier-Stokes 
equations. 9 

In the past, the modal equations were used to model structures for the purpose of 
aeroelastic analysis. For simple geometries such as clean wings, the modal approach can 
produce accurate response results. However, the modal approach may be less accurate 
for complex structures such as wing-body configurations. In order to accurately represent 
aeroelastic responses of general wing-body configurations, the modal equations should be 
replaced with the finite element equations. Recently, a typical wing-body configuration 
has been used to demonstrate aeroelastic responses at transonic Mach numbers using 
the Navier-Stokes equations for fluids and the finite element equations for structures. 10 
Simple one- dimensional beam elements axe used to model the wing-body structures. Each 
node has three degrees of freedom (DOF) corresponding to transverse displacement and 
to transverse and torsional rotations, respectively. 

Recently, a version of ENSAERO 11 that uses the Euler equations for fluids and the 
modal equations for structures has been parallelized on the Intel iPSC/860 at Ames. The 
Intel iPSC/860 is a distributed-memory, multiple-instruction, multiple-data (MIMD) com- 
puter with 128 processors. In this parallel implementation, a domain decomposition ap- 
proach is used in which the fluid equations and the structural equations are modeled in 
separate computational domains. Each domain is mapped individually onto a group of 
processors, referred to as a cube on the Intel iPSC/860. However, because of the coupling 
between the disciplines, there is a need to exchange data, such as pressures and structural 
deformations at interfaces. This exchange between the fluid and structural domains is 
accomplished through an intercube communication mechanism, 12 which enables different 
processors in each cube to communicate directly. 
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III. CURRENT WORK 


1. Finite Element Modeling of Structures 


The finite element representation of structures generally provides more accurate mod- 
eling of structures than the modal representation does in aeroelastic computations. For 
the aeroelastic computations of wing and wing-body structures, plate and shell models 
axe used as shown in Appendix A and C. The parallel implementation of the structural 
domain is mainly discussed in Appendix B. Although the current implementation is only 
using ANS4 plate/shell elements to model structures, it is possible to use different types 
of elements. 

2. Fluid-Structural Interface 


In aeroelastic analysis, it is necessary to represent the equivalent aerodynamic loads at 
the structural nodal points and to represent the deformed structural configurations at the 
aerodynamic grid points. In the present domain decomposition approach, coupling between 
the fluid and structural domains is achieved by interfacing the boundary data, such as 
aerodynamic pressures and structural deflections, at each time step. An analytical moving- 
grid technique has been successfully used to deform the aerodynamic grid according to the 
structural deflections at the end of every time-step. 3 ’ 4 ’ 10 There are different approaches 
for obtaining the external load vector, depending on the equations used for the structural 
dynamic analysis. 

In order to replace modal equations with finite element structural dynamics equations 
for fluid-structural interaction problems, a new fluid-structural interface, similar to the 
modal matrix used for modal equations, should be developed. Several numerical proce- 
dures have been developed for exchanging the necessary information between the fluid and 
structural domains. 13-15 In this research, two different types of fluid-structural interfaces 
were studied and compared as shown in Appendix A. 

The current implementation of the fluid-structural interface on the Intel iPSC/860 
is based on direct node-to-node communication. Thus each of precessors assigned to the 
fluid domain can communicate with any processors of the structural domain. Figure 1 
show the typical communication patern obtained during the aeroelastic computation of a 
High Speed Civil Transport (HSCT) model on the Intel iPSC/860. More details are found 
in Appendix C. 

3. Parallel Integration 


In a serial computer, the integration of both fluid and structural equations is performed 
one after the other in a sequential nature. Figure 2(a) shows the sequential integration 
scheme on serial computers. When implementing the integration scheme on parallel com- 
puters, all processors can be used to solve the fluid and structural equations sequentially 
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as shown Fig. 2(b). But this approach requires more memory per processor and two dis- 
ciplines have to be implemented in a single program. As a result, modularity of each 
algorithm for individual disciplines will have to be sacrificed to a significant degree. In 
addition, this approach will be less efficient as increasing the number of processors because 
the problem is not linearly scaled. 

However, while keeping modularity of each discipline, computations can be done more 
efficiently on MIMD parallel computers by executing the integration of both fluid and 
structural equations concurrently as shown in Fig. 2(c). In the proposed parallel integration 
scheme, both domains start computations independently and one of the solvers waits until 
the other finishes its calculation. Then they exchange the required data with each other 
for the next time step. By doing so, the parallel integration can reduce the idle time 
since only one cube (the fastest) will have to wait. This integration scheme exploits 
the parallelism offered by the domain decomposition approach to solve the coupled fluid- 
structural interaction problems. 

4. Computational Results 


Aeroelastic computations for the NASA Langley Clipped Delta Wing were performed 
first on a Cray Y-MP serial computer and then on the Intel iPSC/860 MIMD parallel 
computer. Figure 3 shows highly deformed wing due to the fluid-structural interactions. 
The surface is colored by pressure coefficient. Comparison between the modal and finite 
element analyses for structures and accuracy of the interfaces considered in this research 
are given in Appendix A. 

Computational performance results are found in Appendix B. Including the data ex- 
change between fluid and structural domains, the current aeroelastically deforming grid 
scheme requires about 12 percent of the computational time per each integration step. 

Application of the procedure to wing-body configurations can be found in Appendices 
A and C. Aeroelastic computations are done on flexible wing and body structures. Figure 
4 shows highly deformed Boeing 1807 HSCT model due to the fluid-structural interactions. 
The surface is colored by pressure coefficient. 


IV. SUMMARY AND CONCLUSIONS 

In this work, a procedure to compute aeroelastic responses on MIMD parallel comput- 
ers using direct coupling of the finite difference Euler/Navier Stokes equations for fluids and 
the finite element dynamics equations for structures is developed for wing and wing-body 
structures. The procedure is based on a domain decomposition approach which enables 
algorithms for the fluid and structural disciplines to be developed and maintained inde- 
pendently. This approach provides an efficient and effective environment to researchers. A 
researcher working in the fluid or the structural discipline can develop his own algorithms 
independent of the others. The only thing to be done together is coupling of the disci- 
plines. Since coupling of the disciplines is achieved by exchanging boundary data through 
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an intercube communication mechanism that does not interfere with interprocessor com- 
munication within a cube, coupling should not cause any problem, This makes it easy 
for each discipline to incorporate and develop new algorithms or data structures without 
interferences. 

The performance of the structural domain is far behind that of the fluid domain. This 
is due to the less desirable performance of the JPCG algorithm. It is noted that direct 
solvers are still in the early stages of development. However, since the procedure developed 
here allows for one domain to select algorithms independent of others, the JPCG algorithm 
can be easily replaced with more efficient algorithms when available. Although the solver 
for the structure is not efficient on a serial computer, reasonable computational speed and a 
good load balance can be achieved by assigning more processors to the structural domain. 
The overall time per integration step of parallel ENSAERO using 96 processors on the 
iPSC/860 is reduced to about 60% of the best time obtained on a single Y-MP processor 
for the particular problem considered. This shows the advantage of using the domain 
decomposition approach for the multidisciplinary analysis on MIMD parallel computers. 

In addition, a parallel integration procedure is developed to solve fluid and structural 
equations together. The parallel integration scheme enables the combination of advanced 
CFD and CSD technologies with minimal increase in computational time per integration 
step while keeping modularity of each discipline. The time per integration step is solely 
determined by the domain that requires most computational time on the iPSC/860. This 
parallel integration is one of the advantages of using MIMD computers for multidisciplinary 
analysis. The procedure developed in this research will provide an efficient tool for solving 
aeroelastic problems of complete aerospace vehicle configurations on MIMD computers. 

Based on this work the following conclusions can be made. 

1. It is feasible to directly couple the finite difference flow equations and finite 
element structural equations to obtain accurate results, though each discipline 
is solved in a separate computational domain. This domain decomposition ap- 
proach takes advantage of efficient methods developed for each individual disci- 
pline. 

2. The use of finite element structures in place of modal structures produces more 
accurate and detailed results. 

3. This domain decomposition approach is suitable for parallel computers. However, 
the structural domain requires a more efficient solver for the application of larger 
size problems. 
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Figure 1 . An intercube communication pattern between fluid and structural domains 
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Figure 2. Comparison between integration schemes 
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Figure 3. Aeroelastic deformation of NASA Langley Clipped Delta wing with pressure contour 
distribution. Mach number = 0.854, Angle of attack = 0.0 degree, Euler Computation 





Figure 4. Aeroelastic deformation of Boeing 1807 HSCT model with surface pressure 

distribution. Mach number = 2.1 , Angle of attack = 4.75 degrees, Euler computation 
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Abstract 

A computational procedure is presented to study 
fluid-structural interaction problems for three-dimen- 
sional aerospace structures. The flow is modeled using 
the three-dimensional unsteady Euler/Navier-Stokes 
equations and solved using the finite-difference ap- 
proach. The three dimensional structure is modeled 
using shell/plate finite-element formulation. The two 
disciplines are coupled using a domain decomposition 
approach. Accurate procedures both in time and space 
are developed to combine the solutions from the flow 
equations with those of the structural equations. Time 
accuracy is maintained using aeroelastic configuration- 
adaptive moving grids that are computed every time 
step. The work done by aerodynamic forces due to 
structural deformations is preserved using consistent 
loads. The present procedure is validated by comput- 
ing the aeroelastic response of a wing and comparing 
with experiment. Results are illustrated for a typical 
wing-body configuration. 

Introduction 

In recent years, significant advances have been 
made for single disciplines in both computational fluid 
dynamics (CFD) using finite-difference approaches 1 
and computational structural dynamics (CSD) using 
finite-element methods (see Chapter 1 of Ref. 2). For 
aerospace vehicles, structures are dominated by inter- 
nal discontinuous members such as spars, ribs, panels 
and bulkheads. The finite-element (FE) method, which 
is fundamentally based on discretization, has proven to 
be computationally efficient to solve aerospace struc- 
tures problems. The external aerodynamics of aerospace 
vehicles is dominated by field discontinuities such as 
shock waves and flow separations. Finite-difference 
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(FD) computational methods have proven to be effi- 
cient to solve such problems. 

Problems in aeroelasticity associated with non- 
linear systems have been solved using both uncoupled 
and coupled methods. 3 Uncoupled methods are less ex- 
pensive but are limited to very small perturbations 
with moderate nonlinearity. However, aeroelastic prob- 
lems of aerospace vehicles are often dominated by large 
structural deformations and high flow nonlinearities. 
Fully coupled procedures are required to solve such 
aeroelastic problems accurately. 

In computing aeroelasticity with coupled proce- 
dures, one needs to deal with fluid equations in an 
Eulerian reference system and structural equations in 
a Lagrangian system. Also, the structural system is 
physically much stiffer than the fluid system. As a re- 
sult, the numerical matrices associated with structures 
are orders of magnitude stiffer than those associated 
with fluids. Therefore, it is numerically inefficient or 
even impossible to solve both systems using a mono- 
lithic numerical scheme. 

Guruswamy and Yang 3 presented a numerically 
accurate and efficient approach to solve this problem 
for two-dimensional airfoils by independently model- 
ing fluids using FD-based transonic small perturbation 
(TSP) equations and st ructures using FE equations and 
coupling the solutions only at boundary interfaces be- 
tween fluids and structures. The coupling of solutions 
at boundaries can be done either explicitly or implicitly. 
This domain decomposition approach allows one to take 
full advantage of numerical procedures of individual 
disciplines such as FD for fluids and FE for structures. 
This accurate coupling procedure has been extended to 
three-dimensional problems and incorporated in sev- 
eral advanced aeroelastic codes such as XTRAN3S 4 , 
ATRAN3S 5 and CAP-TSD 6 based on the TSP the- 
ory. It was later demonstrated that the same technique 
can be used by modeling the fluids with Euler/Navier- 
Stokes equations on moving grids. 7,8 The accuracy of 
the coupling is maintained by matching the field grid 
displacements with the structural displacements at the 
surface. This new development is incorporated in the 
computer code ENSAERO. 9 

As an alternate to the domain decomposition 
approach, there have been some attempts to solve 
both fluids and structures in a single computational 
domain. 1011 This single computational domain ap- 
proach is not new to the researchers dealing with fluid- 
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structural interaction problems. In the late 60 s, there 
were several attempts to solve fluid-structural interac- 
tion problems using a single FE computational domain 
(see Chapter 20 of Ref. 12). The main bottleneck arose 
from ill-conditioned matrices associated with two phys- 
ical domains with large variations in stiffness proper- 
ties. As a result, a subdomain approach was devised 
where fluids and structures are solved in separate do- 
mains and solutions are combined through the bound- 
ary conditions similar to the domain-decomposition ap- 
proach explained above. However, there have been re- 
newed attempts to solve both fluids and structures in 
a single computational domain for aeroelastic applica- 
tions. So far, such attempts are limited to simple two- 
dimensional problems and have not proven to be better 
than the domain decomposition approach. Because of 
the lack of comparison with other approaches and de- 
tails about the computational speed, it is difficult to 
estimate ttie scope of these alternate approaches. The 
drop in the convergence rate from the rigid case to the 
flexible case in Ref. 11 indicates the weakness of the 
single domain approach. 

In the domain decomposition approach, to date, 
advanced CFD methods such as those based on the 
Navier-Stokes equations are used to compute aeroe- 
lasticity of simple wings modeled structural equations. 
The modal approach significantly reduces the number 
of structural unknowns to a great extent when com- 
pared to a direct use of FE equations. For simple ge- 
ometries such as isolated wings, the modal approach 
can produce accurate response results. However, it can 
be less accurate for complex structures such as wing- 
body configurations. Since the structural properties of 
the body are considerably different from those of the 
wing, it is difficult to pre-select the modes to accurately 
represent the full configuration. Therefore, it is more 
accurate to directly use FE structural equations. Also, 
by using the FE equations, stresses and other data that 
are required for the design can be directly computed in 
addition to displacement responses. 

In this work, the capabilities of ENSAERO are 
extended to compute the aeroelastic responses of gen- 
eral wing-body configurations using the Euler/Navier- 
Stokes equations for fluids and plate/shell finite-element 
equations for structures. The coupled equations are 
solved using a time-integration method with configura- 
tion-adaptive moving grids. The results are validated 
for wings and demonstrated for typical wing-body con- 
figurations. Typical aeroelastic responses are computed 
at transonic Mach numbers. 

Governing Aerodynamic Equations 

The strong conservation law form of the Navier- 


Stokes equations is used for shock capturing purposes. 
The thin-layer version of the equations in generalized 
coordinates can be written as 13 

3 T Q + fl € £ + ft,F + d c G= Re~ l d c S (1) 

where Q , E , F, G, and S, are flux vectors in generalized 
coordinates. The following transformations are used in 
deriving Eq. (1). 

r — t 

Z = t(x,y,z,t) 

C = C (x,y,z,t). 

It should be emphasized that the thin-layer approxima- 
tion is valid only for high Reynolds number flows, and 
that very large turbulent eddy viscosities invalidate the 
model. 

To solve Eq. (1), ENSAERO has time-accurate 
methods based on both central-difference and upwind 
schemes. 14 In this paper, the central-difference scheme 
based on the implicit approximate factorization algo- 
rithm of Beam and Warming 15 with modifications by 
Pulliam and Chaussee 16 for diagonalization is used. 
This scheme is first order accurate in time. 

For turbulent flow, the coefficient of viscosity 
needed for Eq. (1) is modeled using the Baldwin- Lomax 
algebraic eddy- viscosity model. 17 All viscous compu- 
tations presented in this paper assume fully turbulent 
flow. This approximation is consistent with the high 
Reynolds number assumption. For vortex- dominated 
flow structures of highly swept wings, a modification to 
the original Baldwin-Lomax model is required. For this 
study, the Degani-Schiff modification 18 to the original 
model for treating vortical flows is used. 

Aeroelastic Equations of Motion 

Following the formulation given in Chapter 20 of 
Ref. 12, the FE matrix form of the aeroelastic equa- 
tion of motion is 

[Mm + [OK?} + [A']{ 9 } = {Z} (3) 

where [M], [G], and [K] are the global mass, damping, 
and stiffness matrices, respectively. {Z} is the aerody- 
namic force vector corresponding to the nodal displace- 
ment vector {q}. 

In this work, it is assumed that the wing-body 
configuration can be modeled using plate/shell ele- 
ments. For this purpose, it is further assumed that 
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the structural properties of the body and wing are rep- 
resented by equivalent shell and plate elements. The 
ANS4 shell/plate element is used to represent the struc- 
tural properties of the wing-body configuration. 19 The 
ANS4 element shown in Fig. 1 is a 20 degrees-of- 
freedom (DOF) element which can model both plates 
and shells. It is based on assumed natural strain ap- 
proach. For the wing-body configuration considered in 
this work, the wing and the body are modeled using 
plate and shell options of the ANS4 element, respec- 
tively. At each node, the DOF allowed are the inplane 
displacements u and u, transverse deflection (w), rota- 
tion about x-axis (0) and rotation about y-axis (<j>). 

The main effort after selecting the FE model of 
the structure falls into computing the global force vec- 
tor {Zj of Eq. (4). { Z } is computed by solving the 
Euler/Navier-Stokes Eq. (1) at given time, t. First, the 
pressures are computed at all surface grid points. The 
forces corresponding to the nodal DOF are computed 
using the FE nodal fluid-structural interfaces discussed 
in the next section. 

Fluid-Structural Interfaces 

In aeroelastic analysis, it is necessary to represent 
equivalent aerodynamic loads at the structural nodal 
points and to represent deformed structural configura- 
tions at the aerodynamic grid points. In the present 
domain decomposition approach, coupling between the 
fluid and structural domains is achieved by combining 
the boundary data such as aerodynamic pressures and 
structural deflections at each time step. An analyti- 
cal moving grid technique has been successfully used to 
deform the aerodynamic grid according to structural 
deflections at the end of every time step.' There are 
several different ways to obtain the global force vector 
{Z} of Eq. (3) depending on the equations used for the 
structural dynamic analysis. 

A number of numerical procedures have been de- 
veloped to exchange the necessary information between 
the aerodynamic and structural domains. 7 A bi-linear 
interpolation and a virtual-surface interface are used 
in this study. The bi-linear interpolation is also called 
the lumped load (LL) approach. In this approach, the 
force acting on each element of the structural mesh is 
first calculated, and then the element nodal force vector 
is obtained by distributing the total force. The global 
force vector is obtained by assembling the nodal force 
vectors of each element. In addition, the deformed con- 
figuration of the CFD grid at the surface is obtained 
by linearly interpolating nodal displacements at finite- 
element nodes. This approach does not conserve the 
work done by the aerodynamic forces and needs fine 
grids for both fluids and structures to give accurate re- 


sults. 

An alternate to the above LL approach is an im- 
proved approach based on virtual surface (VS). In this 
approach, a mapping matrix developed by Appa 20 is 
selected to accurately exchange data between the fluid 
and structural interface boundaries. The reason for se- 
lecting Appa’s method is that the mapping matrix is 
general enough to accommodate changes in fluid and 
structural models easily. In addition, this approach 
conserves the work done by aerodynamic forces when 
obtaining the global nodal force vector. This method 
introduces a virtual surface between the CFD surface 
grid and the finite element mesh for the wing. This 
virtual surface is discretized by a number of finite el- 
ements, which are not necessarily the same elements 
used in the structural surface modeling. 

By forcing the deformed virtual surface to pass 
through the given data points of the deformed struc- 
ture, a mapping matrix relating displacements at struc- 
tural and aerodynamic grid points is derived as 

m =&a}(r 1 [K] + tt.} T M- 1 [4>,) T (4) 

where 

[K]: the free-free stiffness of the virtual surface 
displacement mapping from virtual to 
structural grids 

ip a : displacement mapping from virtual to 
aerodynamic grids 

6: penalty parameter 

Then, the displacement vector at the aerodynamic grid, 
{q a } } can be expressed in terms of the displacement 
vector at the structural nodal points, q S} as 

{?a} ={T] {<?,}■ 

From the principle of virtual work, the nodal force vec- 
tor, { Z $ }, can be obtained as 

{Z.} = [T] t {Z a } 

where {Z a } is the force vector at the aerodynamic grids. 
This procedure is illustrated in Fig. 2. 

The aeroelastic equation of motion, Eq. (3), is 
solved by a numerical integration technique based on 
the constant- average- acceleration method. 

Aeroelastic Configuration Adaptive Grids 

One of the major difficulties in using the Euler/ 
Navier-Stokes equations for computational aerodynam- 
ics lies in the area of grid generation. For steady flows, 
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advanced techniques such as blocked zonal grids 1 are 
currently being used. However, grid-generation tech- 
niques for aeroelastic calculations which involve mov- 
ing components, are still in the early stages of develop- 
ment. In Ref. 7, aeroelastic configuration adaptive dy- 
namic grids were successfully used for computing time- 
accurate aeroelastic responses of wings using a C-H grid 
topology. 

In this work, an H-0 type grid topology is used 
(H in the streamwise and O in the spanwise directions) 
for wing-body configuration. This type of grid topology 
is more suitable for a genera! wing-body configurations. 
It gives better surface grid resolution on the body when 
compared to the C-H grid topology used in Ref. 7. The 
base surface grid is generated using the S3D code. 21 
From the surface grid, the field grid is generated us- 
ing an analytical approach. In this approach, grid lines 
in the radial direction away from the surface are gen- 
erated line by line in the planes normal to the x-axis. 
The new grid lines are generated in such a way that 
the radial lines are approximately normal to the previ- 
ous line. For example, the first line from the surface is 
generated such that the radial lines are approximately 
normal to the surface. In this process the spacing be- 
tween lines are exponentially increased away from the 
surface. This base grid is used for computing pressures 
on the rigid configuration. For aeroelastic analysis, the 
displacements at structural nodes are computed first 
using Eq. (3). These displacements are then mapped 
onto the surface grid points by the interface approach 
discussed above. Finally the field grid is analytically 
generated starting from the new deformed surface. 

Results 

Computations on Wing Configuration 

To demonstrate aeroelastic computations, a typ- 
ical fighter type wing of aspect ratio three and taper 
ratio 1/7 with the NACA 65A006 airfoil section was 
selected. The sweep angle at the quarter chord line 
(A c / 4 ) is 45 deg. The transonic flutter characteristics 
of this wing are available from wind tunnel tests 22 for 
various flow parameters. 

In this computation, the flow field is discretized 
using a C-H grid topology of size 151 x 30 x 35. The 
20 DOF ANS4 shell/plate element 19 was used for the 
FE modeling of the wing structure. The wing is mod- 
eled as a flat plate. Considering the wing structure 
used in the experiment, variation of mass density is al- 
lowed along both chordwise and spanwise directions. 
However, the thickness of the finite element model is 
kept constant. This is based on assumptions that the 
stiffness of the wing is dominated by the aluminum- 


alloy insert and the mass distribution of the wing is 
significantly changed due to plastic foams covering the 
aluminum- alloy insert. This finite-element plate model 
predicts natural vibration modes of the wing that com- 
pare well with the experiment. The first three modal 
frequencies computed by using the finite element model 
are 21.8, 78.1, and 126 Hz and corresponding values 
measured in the experiment are 21.6, 79.7, and 121 Hz, 
respectively. 

This is the first time a shell/plate FE model 
has been directly coupled with the Euler/Navier-Stokes 
equations. As a result, the validity of the coupling ap- 
proach will be verified by comparing the FE results with 
those from the previously well-validated modal analy- 
sis. In this calculation, the FE computations were made 
using 36 plate elements and the modal computations 
were made using the first six modes of the wing. Six 
elements each were assigned along the chordwise and 
spanwise directions, respectively. Figure 3 shows the 
identical displacement responses of the leading edge at 
the tip obtained by both FE and modal analyses for 
A/oo = 0.854, p — 0.70 psi and a = 1.0 deg. Dy- 
namic aeroelastic computations were made setting a 
high value for the damping coefficient so that the final 
results would approach to steady state conditions. The 
VS approach was used to calculate nodal forces for both 
FE and modal analysis. Results in Fig. 3 demonstrate 
the validity of the coupling of plate elements with the 
Euler/Navier-Stokes equations. The FE approach gives 
displacements about 0.1 % higher than the modal ap- 
proach. Such results are expected since the modal ap- 
proach yields a structure that is stiffer than the actual 
one, whereas the FE approach represents the actual 
structural stiffness. 

The accuracy of the results can depend on the 
type of interfaces between fluids and structures. In 
the following calculations the simple lumped load and 
the more accurate virtual surface interfaces are com- 
pared to each other and the results are shown in Fig. 4. 
The wing structure was modeled using 100 ANS4 el- 
ements. Ten elements each were assigned along the 
chordwise and spanwise directions, respectively. For a 
given dynamic pressure of 1.0 psi and initial accelera- 
tion of 1.0 x 10 5 inches/sec, the time history of total 
lift on the wing is presented in Fig. 4. The total lift ob- 
tained by integrating the pressure coefficients at CFD 
grid points is also shown in the figure. The total lift us- 
ing CFD grid points is more accurate than those from 
VS and LL methods. Both VS and LL approaches ob- 
tain the total lift by summing the forces at the FE nodal 
points, which was transformed from the pressure coef- 
ficients through interfaces. The VS approach transfers 
pressure data more accurately than the LL approach. 
The LL approach shows that the response around peaks 
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deviates from the CFD solution. For this case the LL 
approach shows favorable agreement with the VS ap- 
proach. 

Aeroelastic responses were also computed for var- 
ious dynamic pressures in order to predict flutter dy- 
namic pressure and compared with the experment. 22 
Figure 5 shows the stable, near neutrally stable, and 
unstable responses of wing tip displacements at the 
leading edge for dynamic pressures of 0.85, 0.80, and 
0.75 psi for = 0.854. The Navier-Stokes equations 
and the virtual surface interface are used to obtain the 
fE nodal force vector. From the responses shown in 
Fig. 5, the interpolated dynamic pressure for the neu- 
trally stable condition is 0.79 psi. It is noted that the 
experimental dynamic pressure measured at the neu- 
trally stable condition was 0.91 psi. Considering the 
lack of experimental pressure data on the wing and the 
error involved in modeling the wing as a plate with con- 
stant thickness, the result is a favorable prediction of 
the flutter dynamic pressure. 

ENSAERO has capability of modeling both the 
Euler and Navier-Stokes equations. It is of interest to 
know the effect of the type of flow equations on aeroe- 
lastic responses. Such studies will lead to the right 
choice of methods. For this purpose computations are 
made at a high-transonic Mach number of 0.970. Fig- 
ure 6 shows the comparison between the steady pres- 
sures obtained from Euler and Navier-Stokes solutions. 
Since the Mach number is high-transonic, viscous ef- 
fects are dominant. As a result there are significant dif- 
ferences between the Euler and Navier-Stokes solutions 
near and behind the shockwave. The Navier-Stokes so- 
lutions predict lower negative pressures near the shock 
waves. The viscous effects on the integrated total lift 
is shown in Fig. 7. The influence of viscous effects on 
the aeroelastic responses are shown in Fig. 8. In this 
case, aeroelastic computations are made when the wing 
is pitching up to one degree angle of attack (AoA) at a 
pitch rate of 0.01. Because of the reduced aerodynamic 
loads, the tip response from the Navier-Stokes solution 
is lower than that from the Euler solution. 

Computations on Wing-Bodv Configuration 

As stated in the introduction one of the main 
reasons for modeling structures directly by finite ele- 
ments instead of modes is to extend the fluid/structure 
interaction computational capability to more complex 
structures. The procedures demonstrated in the pre- 
vious section are not limited to simple wing configu- 
rations. In this section, results are demonstrated for 
general wing-body configurations where it is not trivial 
to pre-select modes. 


The selected wing-body configuration shown in 
Fig. 9 was modeled using a H-0 type grid topology us- 
ing a grid size of 99 x 79 x 30. Earlier work indicated 
that this grid was adequate for transonic flow compu- 
tations at moderate angles of attack. 23 

In order to study the effect of structural flexibil- 
ity on the flow, aeroelastic computations were made for 
the above wing-body configuration. Both the body and 
wing are allowed to be flexible. The wing is modeled 
using 30 plate elements and the body is modeled using 
90 shell elements. The FE layout is shown in Fig. 10. 
Symmetric boundary conditions are applied at the top 
and bottom body symmetry lines. All DOF are con- 
strained along the wing-body junction. This results in 
a total of 646 DOF for structures. This FE capability 
is incorporated in ENSAERO Version 3.1 in a modular 
way. The skyline data structure is used for the global 
stiffness and mass matrices. 

The structural properties required for the analysis 
results in frequencies that represent a typical transport 
type wing-body. Figure 11 shows the mode shapes of 
the first four modes. For the current structural prop- 
erty assumptions, the first four modes are dominated 
by wing modes. The present 148-node FE model of the 
wing-body configuration can compute up to 646 modes. 

As stated earlier, an analytical moving grid ca- 
pability is implemented in ENSAERO based on H-0 
topology. The grid generated by the code when both 
the wing and the body are deformed is shown in Fig. 12. 
It is noted that the singular planes upstream of the 
leading edge and downstream of the trailing edge are 
deformed according to the deformed shape of the con- 
figuration. 

Forced Motion of Flexible Configuration 

In order to verify the coupling of the surface move- 
ment with the grid movement, computations are made 
by forcing the motion. Computations are made at 
- 0.90, a = 0.0 deg and a reduced frequency k(= ujc/U) 
equal to 0.50, allowing the configuration to deform in 
the first torsional mode of the wing. The wing un- 
dergoes a torsional mode such that the maximum tor- 
sional angle at the tip is 1 deg. The unsteady com- 
putations are started from the converged steady state 
solution and 2400 time steps per cycle of oscillation 
are required. This corresponds to a nondimensional 
computational time step size At = 0.0058. Figure 13 
shows the wing sectional lifts for various sections. As 
expected, the magnitude of the sectional lift increases 
towards the tip. A periodic lift response is obtained 
within two cycles of oscillations. 
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Free Motion on Aeroelastic Configuration 

In this section, aeroelastic computations are made 
on the flexible wing-body configuration by directly cou- 
pling the pressures computed solving the Navier-Stokes 
equations with the FE structural equations. The LL 
interface is used for this computation. The structural 
properties of the wing-body configuration are selected 
to represent a typical aircraft. It is assumed that the 
wing-root is 256 inches long and aeroelastic computa- 
tions are made at a dynamic pressure of 1.0 psi. 

Demonstration computations are made for a static 
aeroelastic case when the configuration is ramping up 
from 0 to 5 deg AoA at = 0.90. The configuration 
is pitched up about the axis perpendicular to the wall 
and located at the leading edge of the wing- root. Start- 
ing from the steady state solution the configuration is 
pitched up at a rate of 0.0012 deg per time step. This 
pitch rate was adequate to obtain a stable and accu- 
rate solution. At every time step the static equilibrium 
position is obtained by solving the static aeroelastic 
equations. At the end of each time step a new field 
grid is generated that conforms to the deformed sur- 
face. Figure 14 shows the response of the leading edge 
of the tip section. 

Computational Resources 

The current Navier-Stokes version of ENSAERO 
runs at 380 M FLO PS on the CRAY C90 at Ames Re- 
search Center. To run a rigid case, the code requires 
33 words of central memory per grid point and 7 mi- 
croseconds of CPU time per time step per grid point. 
For the flexible case there is an additional memory re- 
quirement of 1000 words per node and CPU time of 
25 microseconds per time step per node. A typical dy- 
namic aeroelastic response such as that shown in Fig. 14 
requires about 4 CPU hours and 8 million words of cen- 
tral memory. 

Discussions and Conclusions 

A domain-decomposition computational proce- 
dure is developed to compute aeroelastic responses us- 
ing the Navier-Stokes flow solutions directly coupled 
with finite-element structural equations. The proce- 
dure is demonstrated using plate/shell finite elements. 
Aeroelastic computations are made for a typical wing- 
body configuration. Based on this work the following 
conclusions can be made. 

1. It is feasible to directly couple the finite-difference 
flow equations and finite-element structural equa- 
tions to obtain accurate results, though each dis- 
cipline is solved in a separate computational do- 
main. This domain decomposition approach takes 


advantage of efficient methods developed for each 
individual discipline. 

2. The use of finite-element structures in place of 
modal structures produces more accurate and de- 
tailed results. 

3. There is an increase in the requirement of compu- 
tational time (about 8 % ) and memory for FE 
structures compared to the modal structures (for 
modal structures memory requirement is negligi- 
ble). 

4. The present domain decomposition approach will 
be extended for non-linear structures. 

5. This approach is suitable for parallel computers. 
Work is in progress at the Ames Research Center 
to implement ENSAERO on Intel iPSC/860 par- 
allel computer under NASA’s High Performance 
Computing and Communications (HPCC) Pro- 
gram. 
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Fig. 3 Validation of finite element implementation in 
ENSAERO 



Fig. 4 Comparison of lumped load and virtual surface 
interfacing methods. 
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Fig. 5 Aeroelastic responses using a finite element 
model for structures. 
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Fig. 6 Comparison of transonic steady pressures from 
Euler and Navier-stokes solutions. 



Fig. 7 Comparison of total lift from Euler and Navier- 
stokes solutions. 
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Fig. 8 Comparison of dynamic aeroelastic responses 
from Euler and Navier-stokes solutions 
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Fig. 9 Typical wing-body configuration with portions 
of surface and field physical grids. 


Fig. 11 First four natural modes of the wing-body 
configuration. 




Total elements = 120 
Total nodes = 148 
B.C.’s: 1) symmetric conditions 
at top and bottom of 
the body. 

2) fix all DOF s along 
wing-body junction. 

Total eq'ns: 646 


Fig. 10 Finite-element modeling of the wing-body con- 
figuration. 



Fig. 12 Deformed surface and field grids. 
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Fig. 13 Comparison of sectional lift responses for wing 
twisting mode. 



Fig. 14 Static aeroelastic displacements of wing-tip 
leading edge 1 
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Ames Research Center 

SUMMARY 

This paper presents a procedure for computing the aeroelasticity of wings 
on parallel multiple-instruction, multiple-data (MIMD) computers. In this proce- 
dure, fluids are modeled using Euler equations, and structures are modeled using 
modal or finite element equations. The procedure is designed in such a way that 
each discipline can be developed and maintained independently by using a domain 
decomposition approach. In the present parallel procedure, each computational 
domain is scalable. A parallel integration scheme is used to compute aeroelastic re- 
sponses by solving fluid and structural equations concurrently. The computational 
efficiency issues of parallel integration of both fluid and structural equations are 
investigated in detail. This approach, which reduces the total computational time 
by a factor of almost 2, is demonstrated for a typical aeroelastic wing by using 
various numbers of processors on the Intel iPSC/860. 

INTRODUCTION 

In recent years, significant advances have been made for single-discipline use 
of both computational fluid dynamics (CFD), using finite-difference approaches, and 
computational structural dynamics (CSD), using finite-element methods. In single 
disciplines, computations have been made on complete aircraft. However, only a 
limited amount of work has been completed in coupling these two disciplines for 
multidisciplinary applications. The prime reason is the lack of computational power 
for combining the two major computational fields. The development of a new gen- 
eration of parallel computers can possibly alleviate the restriction of computational 
power. 

A multidisciplinary code for computing unsteady flows and aeroelastic re- 
sponses of aerospace vehicles, ENSAERO, has been developed on serial supercom- 
puters at the Computational Aerosciences Branch of the NASA Ames Research 
Center (ref. 1). This multidisciplinary code computes unsteady aerodynamic re- 
sponses of aircraft using the Euler /Navier- Stokes equations. The modal or finite 
element equations are used to model structures. An aeroelastic shape-conforming 
moving grid is used to include the effect of structural deformations on unsteady 
flows. This code is designed in a modular fashion to adopt several different nu- 
merical schemes suitable for accurate aeroelastic computations. The basic coding 
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of ENSAERO can accommodate zonal grid techniques for efficient modeling of full 
aircraft. 

An early version of ENSAERO (ref. 2) has been successfully applied in com- 
puting aeroelastic responses of a rectangular wing by using the Euler equations 
for fluids and the modal equations for structures. The result demonstrates that 
the code can accurately predict the flutter dynamic pressure of a rectangular wing. 
The code was extended to compute aeroelastic responses using the Navier-Stokes 
equations for fluids (ref. 3). Later, it was updated by utilizing an upwind algorithm, 
and the code has been applied to fighter wings undergoing unsteady motions (refs. 4 
and 5) at moderately large angles of attack. This code also has a capability of mod- 
eling moving control surfaces (ref. 6). Furthermore, ENSAERO has demonstrated 
the capability to simulate transonic flows on wing-body configurations using the 
Navier-Stokes equations (ref. 7). 

For simple geometries such as clean wings, the modal approach can produce 
accurate response results. However, the modal approach may be less accurate for 
complex structures such as wing-body configurations. In order to accurately rep- 
resent aeroelastic responses of general wing-body configurations, the finite element 
equations for structures have been implemented in ENSAERO. A typical wing- 
body configuration has been used to demonstrate aeroelastic responses at transonic 
Mach numbers using the Navier-Stokes equations for fluids and the finite element 
equations for structures (ref. 8). Two-noded beam elements are used to model the 
wing-body structures. Each node has three degrees of freedom (DOF) corresponding 
to transverse displacement and to transverse and torsional rotations, respectively. 

In the past, all computations were accomplished serially, on computers such 
as the Cray Y-MP at Ames Research Center. Currently, a version of ENSAERO 
(ref. 9) that uses the Euler equations for fluids and the modal equations for struc- 
tures has been parallelized on the Intel iPSC/860 at Ames. The Intel iPSC/860 is 
a distributed-memory, multiple-instruction, multiple-data (MIMD) computer with 
128 processors. In this parallel implementation, a domain decomposition approach 
is used in which the fluid equations and the structural equations are modeled in sep- 
arate computational domains. Each domain is mapped individually onto a group of 
processors, referred to as a cube on the Intel iPSC/860. As a result, each discipline 
can be developed and implemented independently of the others. However, because 
of the coupling between the disciplines, there is a need to exchange data, such as 
pressures and structural deformations at interfaces. This exchange between the 
fluid and structural domains is accomplished through an intercube communication 
mechanism (ref. 10), which enables different processors in each cube to communicate 
directly. 

In this work, procedures to compute aeroelastic responses on MIMD parallel 
computers using direct coupling of the Euler equations for fluids and the modal or 
finite element equations for structures are investigated for wings. The implementa- 
tion of the structural domain on the iPSC/860 is described in detail. In addition, the 
computational efficiency issues of parallel integration of both fluid and structural 
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equations axe investigated in detail. The proposed integration scheme exploits the 
architecture of MIMD computers. This research will provide an efficient procedure 

for aeroelastic analysis on MIMD computers. 

This work was completed using the resources of the Numerical Aerodynamic 
Simulation (NAS) Program at NASA Ames Research Center. The work by C. Byun 
was supported by NASA Ames Research Center under Cooperative Agreement 
Number NCC2-740. 

GOVERNING AERODYNAMIC EQUATIONS 

The strong conservation-law form of the Euler equations is used for shock- 
capturing purposes. The Euler equations in generalized coordinates can be written 
as (ref. 11) 

d r Q + d i E + d v F + d < G = 0 (1) 

where Q, E, F, and G are flux vectors in generalized coordinates. The following 
transformations are used in deriving equation (1): 


T = t 

^ = ^(x, j/, z, t) ^) 

T] = rj(x,y,z,t) 

C = C (®» 2/ > t ) 

To solve equation (1), ENSAERO has time-accurate methods based on both 
central-difference and upwind schemes (ref. 12). In this work, the central-difference 
scheme based on the implicit approximate factorization algorithm of Beam and 
Warming (ref. 13) with modifications by Pulliam and Chaussee (ref. 14) for diago- 
nalization was used. This scheme is first order accurate in time. 

AEROELASTIC EQUATIONS OF MOTION 

The governing aeroelastic equations of motion for structures can be written 

(MH$} + (CM + MW = {Z} < 3 > 

where [M], [C], and [K] are the global mass, damping, and stiffness matrices, 
respectively, and where {Z} is the aerodynamic force vector corresponding to the 
displacement vector {q}. These quantities can be expressed in modal coordinates 
or finite element coordinates depending on the method used to obtain structural 
dynamic responses. One of the main efforts is concerned with the computation of 
the global force vector {Z} of equation (3). In this work, for a given time t, {Z} 
is computed by solving the Euler equations. From the solution of equation (1), 
pressure coefficients are computed at all grid points on the wing surface. Using 
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these pressure coefficients, the force vector {Z} is calculated by means of the modal 
matrix or a fluid-structural interface, which is described in the next section. 

The aeroelastic equations of motion (3) have been solved in past work (refs. 2 
and 8) by a numerical integration technique based on the linear acceleration method 
(ref. 15). This method has been successfully used for integrating the modal equa- 
tions, assuming a linear variation of the acceleration. However, this integration 
method requires a very small time step in order to integrate the finite element 
equations of motion. As a result, for the finite element equations of motion, the 
constant-average-acceleration method (ref. 16) is adopted to increase the time-step 
size. This method is an extension of the linear acceleration method. Assuming a 
constant average acceleration on a time interval, the velocities and displacements 
are obtained at a time t as 


{?}t — {?}t-At + 
{?}t = + 


At .... At .... 

~2 W/t-At + ~ 2 ~ 


At {<j}t-At + 


At 2 


({4f}t + At) 


(4) 


Using these equations, the displacements at the end of a time interval can be ob- 
tained by solving 

[D]{q} t = {Z} t + [M]{a}+ [C]{ v} (5) 

where 

4 4 

W = ^2 {?}t-At + ^{<?}t-At + {$’}t-At 

2 

{ v } = ^{tfh-At + 

This is an unconditionally stable scheme, whereas the linear acceleration method is 
conditionally stable. 




FLUID-STRUCTURAL INTERFACES 


In aeroelastic analysis, it is necessary to represent the equivalent aerody- 
namic loads at the structural nodal points and to represent the deformed structural 
configurations at the aerodynamic grid points. In the present domain decomposition 
approach, coupling between the fluid and structural domains is achieved by inter- 
facing the boundary data, such as aerodynamic pressures and structural deflections, 
at each time step. An analytical moving-grid technique has been successfully used 
to deform the aerodynamic grid according to the structural deflections at the end 
of every time step (refs. 1, 2, and 8) There are different approaches for obtaining 
the global force vector {Z} of equation (3), depending on the equations used for 
the structural dynamic analysis. 
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For the modal equations of motion, the global force vector can be easily 
obtained in terms of the preselected mode shapes (modal matrix) as 

{Z} = \l>Vl,m T {A){ZC p } (6) 

where [$] is the modal matrix and [A] the diagonal area matrix of the aerodynamic 
control points. The unsteady differential pressure coefficients on the wing surface 
are defined as {AC P }. It is noted that the modal matrix is also used to represent 
the deformation of the wing. 

Solution of the equations of motion based on the finite element discretization 
requires a fluid-structural interface similar to the modal matrix. Several numerical 
procedures have been developed for exchanging the necessary information between 
the fluid and structural domains (refs. 17-19). However, in this study, a linear 
interpolation scheme is first developed for the interface so that coupling of fluid 
and structural equations could be simple for implementation on the new parallel 
computers. This scheme is called the lumped load interface. In this method, the 
force acting on each element of the structural mesh is first calculated and then 
the element nodal force vector is obtained by distributing the force. The global 
force vector is obtained by assembling the nodal force vectors of each element. In 
addition, the deformed configuration of the CFD grid at the surface is obtained by 
linearly interpolating nodal displacements at finite element nodes. 

Next, a mapping matrix developed by Appa (ref. 18) is selected to accurately 
exchange data between the fluid and structural interface boundaries. The reason 
for selecting Appa’s method is that the mapping matrix is general enough to accom- 
modate changes in fluid and structural models easily. In addition, this approach 
conserves the work done by aerodynamic forces when obtaining the global nodal 
force vector. This method introduces a virtual surface between the CFD surface 
grid and the finite element mesh for the wing. The virtual surface is discretized by 
a number of finite elements, which are not necessarily the same elements used in 
the structural surface modeling. This method is called the virtual surface interface. 

By forcing the deformed virtual surface to pass through the given data points 
of the deformed structure, a mapping matrix relating displacements at structural 
and aerodynamic grid points is derived as 

[T] =[if a ](s~ 1 [^ + M T M' 1 M T W 

where 

[K] is the free-free stiffness of the virtual surface 
ip 3 is the displacement mapping from virtual to structural grids 
ip a is the displacement mapping from virtual to aerodynamic grids 
6 is the penalty parameter 
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Then, the displacement vector at the aerodynamic grid {g 0 } can be expressed in 
terms of the displacement vector at the structural nodal points q a as 

{Qa} = [T] {U 

From the principle of virtual work, the nodal force vector {Z s } can be obtained as 

{ z '.} = m T {z„} 

where { Z a } is the force vector at the aerodynamic grids (ref. 20). 

PARALLEL IMPLEMENTATION OF THE 
AEROELASTIC EQUATIONS 

The domain decomposition approach used in this study enables data struc- 
tures and solution methods for fluid and structural equations to be developed inde- 
pendently. Fluid and structural equations are modeled in separate computational 
domains. However, coupling of the disciplines requires the exchange of the interface 
boundary data, which is accomplished through an intercube communication mech- 
anism (ref. 10). This intercube communication facility enables different processors 
in each cube to communicate directly on the iPSC/860. It is important to keep the 
specification of data exchange routine the same on both computational domains. 

The domain for fluids is capable of solving the Euler equations using 3-D 
uni-partitioning of the computational domain. The uni-partitioning scheme denotes 
that one grid subdomain is assigned to each of the processors. The arrangement 
of processors is described in figure 1. The arrows denote bi-directional data com- 
munication. There are a variety of concurrent algorithms available for solving the 
system of equations for fluids. Currently, the solver for the fluid equations can 
use three different concurrent algorithms: (1) complete exchange-based implemen- 
tations (CE-GE), (2) pipelined-Gaussian elimination (PGE), and (3) substructured 
Gaussian elimination followed by solution of the reduced system by means of bal- 
anced odd-even cyclic reduction (SGE-BCR) (ref. 21). In this work, the one-way 
PGE scheme is used. The choice of algorithms was made largely on the basis of 
memory use. The one-way PGE method allows the use of larger computational 
grids, or of fewer processors, than do the other schemes. More details about the 
implementation of the fluid domain can be found in reference 21. 

For the structural domain, modal equations were first used on the iPSC/860. 
Since a limited number of preselected mode shapes were used, only a single processor 
was assigned to the structural domain. However, in replacing modal equations with 
the finite element equations, it is necessary to use a cube of multiple processors for 
structures. In this study, it is assumed that each subdomain of the entire structure 
is mapped onto a single processor. Each processor stores only the information 
relevant to the subdomain assigned to it. The information can be the stiffness and 
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mass matrices and the applied nodal force vector of a subdomain. Then equation (5) 
is expressed at the subdomain level. 

In this work, a regular finite element mesh is used to model wings as a plate, 
and the domain decomposition is made by using 2-D uni-partitioning as shown 
in figure 1. This type of domain decomposition enables an efficient and simple 
message communication mechanism within the structural domain. Only chordwise 
and spanwise bi-directional messages are exchanged along the subdomain interfaces 
in the structural domain. 

At present, the solver for the structural domain is based on a Jacobi- 
preconditioned conjugate gradient (JPCG) algorithm on the Intel iPSC/860. The 
present JPCG algorithm is based on a parallel conjugate gradient algorithm pro- 
posed by Law (ref. 22). The algorithm is described in the appendix for completeness. 
The advantage of Law’s algorithm is that it does not form the global system matri- 
ces. In this method, the multiplication of a matrix by a trial vector, which is the 
major operation of the conjugate gradient algorithm, is performed at the subdomain 
level. The interprocessor communication is confined to the solution phase. 

INTEGRATION SCHEMES FOR COUPLED DOMAINS 

In a serial computer, the integrations of both fluid and structural equations 
are performed one after the other in a sequential nature. Figure 2(a) shows the 
sequential integration scheme implemented on MIMD computers. In the sequential 
integration scheme, the fluid domain has to wait to proceed to the next time step 
until it receives information about structural deformations. The structural domain 
also has to wait for surface pressure data, so both cubes have their own idle times 
while they wait for data communications. However, since the size of the structural 
equations was small for modal analysis, the effect of the waiting time was negligible 
relative to the computational time per integration step. 

On the contrary, if a large number of modes or direct finite element equations 
are used in order to accurately predict dynamic responses of complex structures, 
the computational time per integration step may be increased rapidly. This is due 
to the increase in the idle time when a sequential integration scheme is used on 
the iPSC/860. However, this situation can be avoided by executing the integration 
of both fluid and structural equations concurrently as shown in figure 2(b). In the 
proposed parallel integration scheme, both solvers start computations independently 
and one of the solvers waits until the other finishes its calculation. Then they 
exchange the required data with each other for the next time step. By doing so, the 
parallel integration can reduce the idle time since only one cube may have partial 
idle time. The resulting speedup achieved by the parallel integration scheme is 
theoretically by a factor of almost 2, provided that computational times required 
for the fluid and structural domains are well balanced. 
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COMPUTATIONAL RESULTS 


To demonstrate an aeroelastic computation, a clipped delta wing of aspect 
ratio 3 and taper ratio 1/7 with the NACA 65A006 airfoil section was selected. 
The sweep angle at the quarter chord line (yl c / 4 ) is 45°. The transonic flutter 
characteristics of this wing are available from wind tunnel tests (ref. 23) for various 
flow parameters. 

In this computation, the flow field is discretized by using a C-H grid topology 
of size 151 x 30 x 25. The CFD grid at the root and the upper surface of the wing 
is shown in figure 3. The CFD grid is assigned to 32 processors on the iPSC/860. 
The processors are arranged as a three-dimensional mesh of eight processors in the 
chordwise direction and two processors in each of the spanwise and surface-normal 
directions. This same arrangement for fluids is kept throughout the computations 
so that the performance of the structural domain can be studied in detail. 

A 20-DOF ANS4 shell element (ref. 24) was used for the finite element mod- 
eling of the wing structure. The wing is modeled as a plate. Considering the wing 
structure used in the experiment, variation of mass density is allowed along the 
chordwise and spanwise directions. But the thickness of the finite element model 
is kept constant. This is based on assumptions that the stiffness of the wing is 
dominated by the aluminum-alloy insert and that mass distribution of the wing is 
significantly changed as a result of plastic foams covering the aluminum-alloy insert. 
This finite element plate model predicts natural vibration modes of the wing that 
compare well with the experiment. The first three modal frequencies computed by 
using the finite element model are 21.8, 78.1, and 126 Hz, and corresponding values 
measured in the experiment are 21.6, 79.7, and 121 Hz, respectively. 

Modal Analysis 

The first parallel version of ENSAERO was capable of using the Euler equa- 
tions for fluids and the modal equations for structures. Using this version of 
ENSAERO, aeroelastic responses were computed. In figure 4, the computed gen- 
eralized displacements of the first three modes for the wing are presented. The 
results were obtained for a freestream Mach number (M^) of 0.854 and a given 
dynamic pressure (P) of 0.7 psi. In this calculation, the first six mode shapes of 
the wing were used to predict the structural dynamic responses. Identical results 
were obtained from serial and parallel computers. At this point, it is verified that 
the fluid and structural domains of the parallel ENSAERO and of the intercube 
communication mechanism are properly working. It should be noted that only one 
processor is assigned for the structural analysis since only six mode shapes are used 
to represent the structural properties of the wing. The wall-clock times required 
for each integration step on a single processor of the Cray Y-MP and on the Intel 
iPSC/860 are 1.36 and 3.03 seconds, respectively. 
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The proposed parallel integration scheme is compared with the sequential 
integration scheme used on serial computers. Using both sequential and paral- 
lel integration schemes, aeroelastic responses of the wing were computed on the 
iPSC/860. Aeroelastic responses were computed by simulating experimental con- 
ditions for a freestream Mach number of 0.977 and for a given dynamic pressure 
of 0.65 psi. Results from sequential and parallel integration schemes agree well, as 
shown in figure 5. Since only six modes were used for the modal analysis, the reduc- 
tion in computational time was marginal (less than 2% of the time per integration 
step used in the sequential integration scheme). By increasing the number of modes 
to 50, the sequential integration scheme required 5% more computational time per 
integration step whereas the time for the parallel integration scheme remained the 
same. This trend is more evident as the number of equations increases. 

Finite Element Analysis 

For simple geometries such as rectangular wings, the modal analysis can 
accuratly predict dynamic responses. However, the modal analysis with a limited 
number of preselected mode shapes may be less accurate for complex structures 
such as wing-body configurations. In order to accurately represent aeroelastic re- 
sponses of general aircraft configurations, the modal analysis was replaced with the 
direct finite element analysis in ENSAERO. The finite element equations were first 
tested on the Cray Y-MP version of the code and then were parallelized on the 
Intel iPSC/860. 

The lumped load and virtual surface interfaces on the Y-MP are shown in 
figures 6 and 7. The wing structure was modeled using 100 ANS4 elements. Ten 
elements each were assigned along the chordwise and spanwise directions, respec- 
tively. The time history of total lift on the wing for a given dynamic pressure of 
1.0 psi and initial acceleration of 1.0 x 10 5 inches/sec 2 , is presented in figure 6. 
The exact solution is the total lift obtained by integrating pressure coefficients at 
CFD grid points. Both virtual surface (VS) and lumped load (LL) interfaces obtain 
the total lift by summing the forces at the finite element (FE) nodal points, which 
were transformed from pressure coefficients through interfaces. The virtual surface 
interface transfers pressure data more accurately than the lumped load interface. 
The lumped load interface shows a favorable result although the response around 
peaks deviates from the exact solution. In addition, the tip displacements of the 
wing at the leading edge are presented in figure 7 . The lumped load approach shows 
favorable agreement with the virtual surface approach. 

The lumped load approach was first used as the fluid-structural interface for 
the finite element equations on the iPSC/860. The choice of interfaces was made 
largely on the basis of memory use. The size of the mapping matrix for the virtual 
surface interface becomes too large to fit on a single processor on the iPSC/860 when 
the number of fluid grid points or structural nodal points on the wing increases. 
However, the lumped load approach requires only a small amount of memory to 
identify the location of fluid grid points on a finite element discretization. 
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In order to check the parallel implementation of the finite element equations, 
the tip displacements of the wing were computed, without aerodynamic forces, on 
the Y-MP and the iPSC/860. This wing was modeled using 64 ANS4 elements. For 
this computation, processors were assigned as a 2-D mesh of two processors in the 
chordwise and spanwise directions, respectively, on the iPSC/860. Identical results 
were obtained on both the Y-MP and the iPSC/860. 

For the structural model using finite element equations, the aeroelastic re- 
sponses were obtained using both sequential and parallel integration schemes on 
the iPSC/860. The results are presented in figure 8. The responses were ob- 
tained for a given dynamic pressure of 1.0 psi. The same finite element mesh 
and processor arrangement used in the previous case are used for this computation. 
The two results agree well. The wall-clock times per integration step achieved are 
3.45 and 3.00 seconds by using sequential and parallel integration schemes, respec- 
tively. The speedup is still marginal since the total number of equations is relatively 
small (360 DOF) for the structural dynamic analysis. However, for 256 finite el- 
ements (1360 DOF) with four processors on the structural domain, the wall-clock 
times per integration step are 6.22 and 3.33 seconds by using sequential and par- 
allel integration schemes, respectively. The speedup achieved is 1.87 by using the 
parallel integration scheme. When the computational time between the fluid and 
struct ural domains is balanced, maximum speedup can be achieved. 

The parallel integration scheme enables the combination of advanced CFD 
and CSD technologies with minimal increase in the computational time per integra- 
tion step. The required computational time per integration step is determined by 
both the fluid and structural domains on serial computers. However, using the par- 
allel integration scheme on MIMD computers, the time is solely determined by the 
computational domain that requires more time per integration step. This parallel 
integration is one of the advantages of using MIMD computers for multidisciplinary 
analysis. 

Aeroelastic responses were also computed for various dynamic pressures in 
order to predict the flutter dynamic pressure. Figure 9 shows the stable, near 
neutrally stable, and unstable responses of wing tip displacements at the leading 
edge for dynamic pressures of 0.80, 0.85, and 0.90 psi, respectively. From the 
responses shown in figure 9, the interpolated dynamic pressure for the neutrally 
stable condition is 0.84 psi. It is noted that the experimental dynamic pressure 
measured at the neutrally stable condition was 0.91 psi. Considering the lack of 
experimental pressure data on the wing and the error involved in modeling the 
wing as a plate with constant thickness, the computational result is an acceptable 
prediction of the flutter dynamic pressure. 

Performance 

In order to support multidisciplinary analysis with practical computational 
turnaround times for design work, the computational domain on parallel machines 
must be scalable. This section describes several aspects of performance, including 
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single processor computational rates, domain decomposition strategy, and scalabil- 
ity of the structural domain in ENSAERO on the iPSC/860. It is noted that only 
a linear model is used for the dynamic analysis of the structural domain. The per- 
formance of the fluid domain can be found in reference 21. All performance data 
reported are for 64-bit arithmetic. 

In order to measure the performance of the structural domain on the Intel 
iPSC/860, the identical code, except message passing routines, was run on the 
Y-MP and the averaged time per integration step was obtained. All FLOP rates 
quoted are calculated by comparing the time per integration step on the iPSC/860 
with that on the Y-MP using a single processor. Operation counts from the Cray 
Hardware Performance Monitor are used. 

A single processor of the iPSC/860 was able to house 256 ANS4 elements. 
The Y-MP equivalent MFLOPS obtained is about 4.2 MFLOPS for this size of 
problem, and the corresponding rate is about 77 MFLOPS on a single Y-MP pro- 
cessor. This rate is about 7% of the peak performance of a single processor on the 
iPSC/860. Similar performance was reported by Ryan and Weeratunga for the fluid 
domain (ref. 21). 

The performance of the structural domain in parallel ENSAERO has been 
measured over a wide range of processor numbers and problem sizes as shown in 
figure 10. The speedup relative to the Y-MP is defined as 


speedup = 


tCray 
t Intel 


where tcray and t Intel are the computational time per integration step measured 
on the Y-MP and the iPSC/860, respectively. Only a single processor is used 
to measure tcray on the Y-MP. The open and filled symbols denote the domain 
decomposition which results in the minimum and maximum bandwidths of the 
stiffness matrix of each subdomain for a given number of processors. 

For the case of 1,360 DOF, the computational time per integration step on 
the iPSC/860 is barely closed to that on the Y-MP when 64 processors are in use. 
However, as increasing the size of problem (10,560 and 20,800 DOF), the iPSC/860 
achieves about the speed of the Y-MP by using 16 processors. It is evident that 
the JPCG solver on the iPSC/860 performs better as the size of problem increases. 
For the case of 20,800 DOF, the relative speedup achieved is about 8 by the time 
64 processors are in use. 

For a given number of processor, the obtained speedup varies depending on 
the domain decomposition strategy as shown in figure 10. Only the results for the 
minimum and maximum bandwidths are presented for clarity. The speedup in- 
creases as decreasing the matrix bandwidth of each subdomain for a given number 
of processors on the iPSC/860. This is due to the fact that the conjugate gradient 
algorithm is subjected to the multiplication of the coefficient matrix and a trial vec- 
tor. Since the multiplication is performed only at the subdomain level in the JPCG 
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solver, the smaller bandwidth results in fewer operations and quicker computational 
time. 

The overall performance of ENSAERO on both the Y-MP and the iPSC/860 
is shown in figure 11 for the case of 113,250 grid points for the fluids and 10,560 DOF 
for the structure. In this computation, 32 processors are assigned to the fluid domain 
and 16 to 64 processors to the structural domain. Both the skyline reduction and 
JPCG solvers are compared on the Y-MP but only the JPCG solver is used for the 
structural domain on the iPSC/860 at the present time. A parallel version of the 
skyline reduction solver is under implementation. 

The height of each column is the time per integration step. Each column is 
divided into zones representing the time spent for the fluid domain, the structural 
domain, and for idle/intercube communication. It should be noted that the time per 
integration step for the skyline reduction solver included only time spent for forward 
reduction and back substitution without the factorization time. The reason is that 
the contribution of the factorization time to the computation time per integration 
step is negligible when a large number of time steps are required to obtain linear 
dynamic responses. Providing that 7,000 steps are required for a typical aeroelastic 
computation of the given wing, the increase of the time per integration step is about 
0.5% of the time used for structure. 

It is evident that the skyline reduction solver outperforms the JPCG solver 
on the Y-MP. However, the JPCG solver is first implemented on the iPSC/860. The 
reason is that it is desirable to compare the performance of the two solvers on the 
iPSC/860. The JPCG solver on the iPSC/860 could not achieve the performance of 
the skyline reduction solver without including the factorization time on the Y-MP. 
However, as far as the JPCG solver is concerned, 16 processors on the iPSC/860 can 
obtain the performance of a single processor on the Y-MP. In addition, the overall 
performance of ENSAERO using 96 processors on the iPSC/860 is about one-t hir d 
of that obtained using the skyline reduction solver with a single processor on the 
Y-MP. This result is based on the averaged time per integration step. It should be 
noted that the structural domain determined the time per integration step for this 
particular problem on the iPSC/860. Most of the time on the fluid domain was spent 
waiting for the interface boundary data. This means that fewer processors can be 
assigned to the fluid domain without sacrificing computational performance as long 
as the memory on each processor can accommodate the assigned grid partitioning. 

CONCLUSIONS 

A parallel version of a multidisciplinary code, ENSAERO, was developed 
on the Intel iPSC/860. A domain decomposition approach was used to enable 
the fluid and structural domains to be developed and maintained independently. 
This approach provides an efficient and effective environment to researchers. A 
researcher concerned with the fluid or the structural domain can develop his own 
discipline independent of the others. The only thing to be done together is coupling 
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of the disciplines. Since coupling of the disciplines is achieved by an intercube 
communication mechanism, coupling should not cause any problem as long as each 
domain maintains the specification for intercube communication. This makes it easy 
for each discipline to incorporate and develop new algorithms or data structures 
without interferences. 

In addition to the modal analysis, the capability of finite element analysis 
for structural dynamic analysis has been added to parallel ENSAERO. The struc- 
tural domain on the iPSC/860 is again divided into a number of subdomains. The 
solution of the structural domain is obtained by the Jacobi preconditioned conju- 
gate gradient (JPCG) solver. The partition of the structural domain for a wing 
is two-dimensional uni-partitioning since the wing is modeled as a plate and the 
discretization is a regular mesh. This enables the message communication within 
the structural domain to be very simple and efficient. 

As far as the structural domain is concerned, the performance on the 
iPSC/860 is still far behind the best performance of the Y-MP a result of the 
poor performance of the JPCG algorithm. A parallel version of the skyline reduc- 
tion solver is being implemented for the structural domain on the iPSC/860. This 
will provide increased performance for the structural domain on the iPSC/860. For 
a problem size of 10,560 DOF, the overall performance of parallel ENSAERO using 
96 processors on the iPSC/860 is about one-third of that obtained using the skyline 
reduction solver for the structural domain on a single Y-MP processor. 

The parallel integration scheme enables the combination of advanced CFD 
and CSD technologies with minimal increase in computational time per integration 
step. The computational time per integration step is solely determined by the do- 
main that requires more computational time on the iPSC/860 whereas that time it 
is determined by both domains on serial computers. This parallel integration is one 
of the advantages of using MIMD computers for multidisciplinary analysis. The pro- 
cedure developed in this research will provide an efficient tool for solving aeroelastic 
problems of complete aerospace vehicle configurations on MIMD computers. 
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APPENDIX 


1. Set diagonal elements of global stiffness matrix 

Send {Kjj} to neighboring processors 
Receive {Kj from neighboring processors 
{<*'} = E {Kjj} + {Kjj} 

2. Set initial trial and residual vectors 

{?'} = 0 
{’■'} = {*•} 

3. Set initial search direction 

Send {r e } to neighboring processors 
Receive {r e } from neighboring processors 
{ 5 '} = E{r‘} + {r'} 

2 j = s'/rfj, j = 1, neq* 
p‘ = 

7o = E p € , e = 1 , np (global sum) 

7 = 7o 
{p e } = {s c } 

4. Operations at subdomain level 

{«■*} = \K'\{p') 

0 ' = {j»*} r {ti*} 

<r e = ff’/y 

5. Update solution and residual 

1/a = £ < 7 C , e = 1 , np (global sum) 

{? e } = {$ e } + <*{p e } 

{r c } = {r c } — a{u e } 

6. Update search direction and check convergence 

Send {r c } to neighboring processors 
Receive {r e } from neighboring processors 
{«*} = £{r<} + {r e } 
z j = s]/d e j, j = l,...,neq e 
P e = {r e } T {z e } 

Inew = £ p e , e = 1, ...,np (global sum) 
IF( inewho < TOLERANCE) STOP 

{p e } = {S e } + (jnew/l) {P e ) 

*1 ~ 'l new 

7. Repeat 4 to 6 until converged 
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Figure 1. Processor arrangements and message exchanges through interprocessor 
and intercube communications. 
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(a) Sequential Integration 



(b) Parallel Integration 


Figure 2. Flow diagrams for sequential and parallel integration schemes. 
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Figure 3. The CFD grid at the root and surface of a clipped delta wing with the 
NACA 65A006 airfoil section. (Aspect ratio = 3.0, Taper ratio = 1/7, A c / 4 = 45°) 
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Figure 4. Comparison among the first three generalized displacement histories ob- 
tained by the serial and parallel ENSAERO codes. (Moo = 0.854, P = 0.70 psi) 



Figure 5. Generalized displacement histories obtained from sequential and parallel 
integration schemes. (Moo = 0.977, P = 0.65 psi) 
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Figure 6. Comparison of total lift responses. 
(Mqo = 0.854, q = 0 deg., P = 1.0 psi) 



Figure 8. Aeroelastic responses obtained 
by using sequential and parallel integra- 
tion schemes. (Moo = 0.854, a = 0 deg., 
P = 1.0 psi) 



Figure 7. Comparison of wing tip displace- 
ments at the leading edge obtained by us- 
ing virtual surface and lumped load inter- 
faces. (Moo = 0.854, a = 0 deg., 
P = 1.0 psi) 



Figure 9. Aeroelastic responses of a clipped 
delta wing at three different dynamic pres- 
sures using Euler equations for the fluid 
and finite element equations for the struc- 
tural domains. (Moo = 0.854, a = 0 deg.) 
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Figure 10. Computational performance of the structural domain in ENSAERO with 
various problem sizes and domain decompositions on the Intel iPSC/860. 
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Figure 11. Overall computational performance of ENSAERO on the Cray Y-MP 
and the Intel iPSC/860. 


22 




APPENDIX C. 






WING-BODY AEROELASTICITY USING FINITE-DIFFERENCE FLUID/ 
FINITE-ELEMENT STRUCTURAL EQUATIONS ON PARALLEL COMPUTERS 

Chansup Byun* and Guru P. Guruswamy** 

Computational Aerosciences Branch 
NASA Ames Research Center, Moffett Field, California 94035-1000 


Abstract 

This paper presents a procedure for computing the 
aeroelasticity of wing-body configurations on multiple- 
instruction, multiple-data (MIMD) parallel computers. 
In this procedure, fluids are modeled using Euler equa- 
tions discretized by a finite difference method, and 
structures are modeled using finite element equations. 
The procedure is designed in such a way that each disci- 
pline can be developed and maintained independently 
by using a domain decomposition approach. A par- 
allel integration scheme is used to compute aeroelas- 
tic responses by solving the coupled fluid and struc- 
tural equations concurrently while keeping modularity 
of each discipline. The present procedure is validated 
by computing the aeroelastic response of a wing and 
comparing with experiment. Aeroelastic computations 
are illustrated for a High Speed Civil Transport type 
wing-body configuration. 

Introduction 

The analysis of aeroelasticity involves solving fluid 
and structural equations together. Both uncoupled 
and coupled methods can be used to solve problems 
in aeroelasticity associated with nonlinear systems. 1 
Uncoupled methods are less expensive but are limited 
to very small perturbations with moderate nonlinear- 
ity. However, aeroelastic problems of aerospace vehicles 
are often dominated by large structural deformations 
and high flow nonlinearities. Fully coupled procedures 
are required to solve such aeroelasticity problems accu- 
rately. 

Such coupling procedures result in an increased 
level of complication. Therefore, aeroelastic analysis 
has been mostly performed by coupling advanced com- 
putational fluid dynamics (CFD) methods with simple 
structural modal equations or advanced computational 
structural dynamics (CSD) methods with simple flow 
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solutions. However, these approaches can be less ac- 
curate for the aeroelastic analysis of practical problems 
such as a full aircraft configuration in transonic regime. 
It is necessary to develop a fully coupled procedure uti- 
lizing advanced computational methods for both disci- 
plines. 

Recently, coupled fluid-structural interaction prob- 
lems are being studied using finite-difference Euler or 
Navier-Stokes flow equations and finite-element struc- 
tural equations of motion as demonstrated by an aeroe- 
lastic code ENSAERO. 2 * 3 However, applications are 
limited to simple structural models. For the compli- 
cated fluid and structural models, computations are 
performed in a step-by-step fashion. 4 The main reason 
is that the use of detailed models for both disciplines 
requires unprecedented computing speeds and amounts 
of memory. The emergence of a new generation of par- 
allel computers can possibly alleviate the restriction on 
the computational power. 

In order to solve the coupled fluid-structural equa- 
tions, some attempts have been made to solve both flu- 
ids and structures in a single computational domain 5,6 . 
The main defect of this approach is the ill-conditioned 
matrices associated with two physical domains with 
large variations in stiffness properties. So far, such 
attempts have been limited to simple two-dimensional 
problems. 

To overcome the difficulties arising from a sin- 
gle domain approach, the domain decomposition ap- 
proach reported in Ref. 1 has been incorporated in sev- 
eral advanced aeroelastic codes such as XTRAN3S, 7 
ATRANS3S 8 and CAP-TSD 9 based on the transonic 
small perturbation theory. This domain decomposition 
approach models fluids and structures independently. 
The coupling of two disciplines is accomplished by ex- 
changing data at interfaces between fluids and struc- 
tures. This allows one to take full advantage of numer- 
ical procedures for individual disciplines such as finite 
difference methods for fluids and finite element meth- 
ods for structures. It was later demonstrated that the 
same technique can be used for modeling the fluids with 
Euler/Navier-Stokes equations on moving grids. 10,11 
The accuracy of the coupling is maintained by match- 
ing the surface grid deformation with the structural 
displacements at the surface. This new development is 
incorporated in the computer code. ENS AERO. 12 Sim- 
ilar work has also been reported recently in Ref. 13. 

For the implementation of the ENSAERO code 



on parallel computers, two types of parallel comput- 
ers were considered. These are the single- instruct ion, 
multiple-data (SIMD) and multiple-instruction, multi- 
ple-data (MIMD) type computers. However, MIMD 
type parallel computers are more suitable for compu- 
tationally efficient implicit solvers and the domain de- 
composition approach used in the ENSAERO code. By 
decomposing the computational domain into a num- 
ber of subdomains and solving an implicit problem on 
each subdomain, a MIMD computer can reduce the in- 
terprocessor communication required for the inversion 
of a large matrix resulting from an implicit method. 
Furthermore, a MIMD parallel computer can exploit 
the parallelism offered by the domain decomposition 
approach for the coupled fluid and structural disci- 
plines; each computational domain can be treated con- 
currently. In addition, each fluid and structural algo- 
rithms can be designed in a modular fashion on MIMD 
parallel computers. 

In this work, a procedure to compute aeroelastic- 
ity on MIMD parallel computers is described. The 
fluid and structural equations on separate computa- 
tional domains are coupled by the exchange of interface 
data. The computational-efficiency issues of parallel in- 
tegration of both fluid and structural equations are in- 
vestigated using a parallel version of ENSAERO. The 
fluid and structural disciplines are modeled using finite- 
difference (FD) and finite-element (FE) approaches, re- 
spectively. The coupled equations are solved using a 
time integration method with configuration-adaptive 
moving grids. The procedure is designed in a modular 
fashion so that each computational discipline can be 
developed independently and be modified easily. The 
aeroelastic computations are demonstrated for a High 
Speed Civil Transport (HSCT) type wing-body config- 
uration on the Intel iPSC/860 parallel computer. 

Aeroelastic Computation 

The governing aeroelastic equations of motion for 
structures can be written as 

[Mm+[cm+[K}{q}= { z } (i) 

where [Af], [C], and [K] are the global mass, damping 
and stiffness matrices, respectively. { Z } is the aerody- 
namic force vector corresponding to displacement vec- 
tor {g}. One of the main efforts is computing the aero- 
dynamic force vector {Z}, which is obtained by solving 
the fluid flow equations. After obt aini ng the aerody- 
namic force, aeroelastic responses can obtained by solv- 
ing eq. (1). A numerical integration technique based on 
the constant-average-acceleration method 14 is used to 
integrate the aeroelastic equations. This is an uncon- 
ditionally stable scheme. 


A domain decomposition approach is selected to 
solve eq. (1) in conjunction with the flow equations. 
Each of the fluid and structural equations is modeled 
in a separate computational domain. Coupling between 
the fluid and structural equations is accomplished by 
exchanging boundary interface data at the end of every 
time step when solving eq. (1). The advantage of this 
approach is that one can select an efficient algorithm 
for the fluid domain regardless of the structural domain 
and vice versa. In this work, a finite difference method 
is selected for fluids and a finite element method for 
structures. 

In the fluid domain the strong conservation law 
form of the Euler equations is used to model the flow. 
To solve the Euler equations, the central-difference 
scheme based on the implicit approximate factorization 
algorithm of Beam and Warming 15 with modifications 
by Pulliam and Chaussee 16 for diagonalization is used. 
The scheme is first order accurate in time. 

To exchange boundary interface data, it is neces- 
sary to represent the equivalent aerodynamic loads (i.e. 
normal stress) at the structural nodal points and to 
represent the deformed structural configurations at the 
aerodynamic grid points. Several numerical procedures 
have been developed to exchange the necessary infor- 
mation between the fluid and structural domains. 17-20 

A node-to-element approach is used to define the 
location of the points of the fluid surface grid relative 
to finite elements at the surface of the structure for 
coupling purposes. In this approach, every grid point 
of the fluid that lies on the fluid-structural interface is 
identified with respect to a finite element as shown in 
Fig. 1. However, in general, it is not straightforward to 
determine the local coordinate information of each grid 
point within a finite element. A numerical inverse map- 
ping technique developed by Murti and Valliappan 21 is 
used to obtain the local coordinate information of all 
the interface points of the fluid grid with respect to 
surface elements of the structure. Once the location of 
each fluid grid point is obtained, the nodal force vec- 
tor can be easily obtained. Also the deformation of the 
fluid surface grid is determined by using shape functions 
of the finite elements used to model the structure. In 
addition, a linear extrapolation is used to compute the 
deformation at the points of the fluid surface grid which 
are not on the surface of the structure, e.g. points at 
singular planes. Starting from the deformed fluid sur- 
face grid, the field grid is generated as explained in a 
later section. 

Parallelization of ENSAERO 

The domain decomposition approach enables data 
structures and solution methods for fluid and struc- 
tural equations to be developed independently. Fluid 



and structural equations are modeled in separate com- 
putational domains. Each domain is mapped individ- 
ually onto a group of processors, referred to as a cube 
on the Intel iPSC/860, which is selected for this work. 
The Intel iPSC/860 is a distributed- memory, multiple- 
instruction, multiple-data (MIMD) computer with 128 
processors. 

Because of coupling between the two disciplines, 
the interface boundary data such as surface pressures 
and structural displacements should be exchanged. 
This exchange between the fluid and structural do- 
mains is accomplished through an intercube commu- 
nication mechanism. 22 This intercube communication 
facility enables different processors in each cube on the 
iPSC/860 to communicate directly. 

The fluid flow algorithm solves the Euler equa- 
tions using 3-D uni-partitioning of the computational 
domain. The uni-partitioning scheme assigns one sub- 
domain grid to each of the processors. The mapping 
of subdomain grids to processors is described in Fig. 2. 
The arrows denote bi-directional data communication. 
There are a variety of concurrent algorithms available 
for solving the system of equations for fluids. More 
details about the implementation of the fluid flow algo- 
rithms can be found in Ref. 23. 

For the structural domain, regular finite element 
meshes are used to model the wing and the body as 
plate and shell structures, respectively. The domain 
decomposition is made by using 2-D uni-partitioning 
as shown in Fig. 2. This type of domain decomposition 
enables an efficient and simple message communication 
mechanism within the structural domain. 

The solver for the structural domain is based on 
a Jacobi-preconditioned conjugate gradient (JPCG) al- 
gorithm on the Intel iPSC/860. The present JPCG al- 
gorithm is obtained by implementing the diagonal pre- 
conditioner to a parallel conjugate gradient algorithm 
proposed by Law. 24 In this method, the structural fi- 
nite element model is divided into subdomains and only 
local matrices related to the subdomains are assembled. 
The multiplication of a matrix by a trial vector, which 
is the major operation of the conjugate gradient algo- 
rithm, is performed at the subdomain level. Interpro- 
cessor communication is confined to the solution phase, 
and the communication is only performed between pro- 
cessors which have common finite element nodes. 

Aeroelastic Configuration Adaptive Grids 

One of the major difficulties in solving the Euler 
equations for computational aerodynamics lies in the 
area of grid generation. For steady flows, advanced 
techniques such as blocked zonal grids 25 are currently 
being used. However, grid-generation techniques for 
aeroelastic calculations, which involve moving compo- 


nents, are still in the early stages of development. Gu- 
ruswamy has developed analytical schemes for aeroelas- 
tic configuration-adaptive dynamic grids and demon- 
strated time-accurate aeroelastic responses of wing 10 
and wing-body 2,3 configurations. 

In this work, an H-0 type grid topology is used (H 
in the streamwise and O in the spanwise directions) for 
wing-body configurations. This type of grid topology is 
more suitable for general wing-body configurations. It 
gives better surface grid resolution on the body when 
compared to the C-H grid topology used in Ref. 10. 
The grid is designed such that flow phenomena such as 
shock-waves, vortices, etc. and their movement around 
the wing-body configurations are accurately simulated. 
The base surface grid is prepared using the S3D code. 26 
From the surface grid, the field grid is generated using 
an analytical approach. In this approach, grid lines in 
the radial direction away from the surface are gener- 
ated line by line in the planes normal to the longitudi- 
nal body axis. First, the radial lines are generated ap- 
proximately normal to the surface. Then, the new grid 
lines in the azimuthal direction are generated in such 
a way that the spacing between lines are exponentially 
increased away from the surface. This method can be 
used for generating the base field grid of the rigid con- 
figuration and the aeroelastically deformed field grid of 
the flexible configuration. 

Parallelization of this approach is accomplished us- 
ing the uni-partitioning scheme in the fluid domain. 
The present approach for aeroelastic configuration- 
adaptive grids only requires the deformed surface grid 
and the coefficients used in the exponential function 
to define the grid spacing between lines away from the 
surface. So the interprocessor co mmuni cation needed 
to generate the deformed field grid within the fluid do- 
main is minimal, and takes place only between proces- 
sors assigned along the surface-normal direction. Each 
of the processors can generate the assigned subdomain 
grid of the deformed field grid concurrently once infor- 
mation about the local surface grid has been broadcast. 

The grid is generated at every time step based on 
the aeroelastically deformed position of the structure. 
First, the displacements at the points of the fluid sur- 
face grid on the structure are obtained on the processors 
assigned to the structural domain. This is done by us- 
ing the local coordinate information and the finite ele- 
ment shape functions. The displacements at the points 
of the fluid surface grid are sent only to the appropriate 
processors on the fluid domain which contain the sur- 
face grid points. Then, a linear extrapolation is used to 
obtain the displacements at the remaining points of the 
surface grid, such as the points at singular planes which 
are not on the surface of the structure. At this stage, 
the deformed surface grid is distributed only to proces- 
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sors of the fluid domain which contain the local surface 
grid points. It should be noted that the deformed sur- 
face grid residing on each processor of the fluid domain 
is only the part of the whole surface grid according to 
the grid partitioning. If two or more processors are as- 
signed along the surface-normal direction, each of the 
deformed surface partitions is sent to processors which 
have the same partitioning indices of the surface grid. 
Finally, all processors of the fluid domain generate their 
subdomain of the deformed field grid concurrently. 

Parallel Integration for Coupled Domains 

In a serial computer, the integration of both 
fluid and structural equations is performed one af- 
ter the other in a sequential nature. Figure 3(a) 
shows the sequential integration scheme implemented 
on MIMD parallel computers. In the sequential inte- 
gration scheme, the fluid domain has to wait to pro- 
ceed to the next time step until it receives information 
about structural deformations. The structural domain 
also has to wait for surface pressure data. So, both 
cubes have their own idle times waiting for data com- 
munications. The computational time per integration 
step will be determined by times spent on both domains 
when a sequential integration scheme is used. In order 
to avoid the idle times between the fluid and structural 
computations, all processors can be used to solve the 
fluid and structural equations sequentially as done in 
serial computations. But this approach requires more 
memory per processor and two disciplines have to be 
implemented in a single program. As a result, modu- 
larity of each algorithm for individual disciplines will 
have to be sacrificed to a significant degree. In addi- 
tion, this approach will be less efficient as increasing 
the number of processors because the problem is not 
linearly scaled. 

However, while keeping modularity of each disci- 
pline, computations can be done more efficiently on 
MIMD parallel computers by executing the integration 
of both fluid and structural equations concurrently as 
shown in Fig. 3(b). In the proposed parallel integra- 
tion scheme, both domains start computations indepen- 
dently and one of the solvers waits until the other fin- 
ishes its calculation. Then they exchange the req uir ed 
data with each other for the next time step. By do- 
ing so, the parallel integration can reduce the idle time 
since only one cube (the fastest) will have to wait. The 
resulting speedup by the parallel integration scheme is 
a factor of almost 2, provided that computational times 
required for the fluid and structural domains are well 
balanced. This integration scheme exploit the paral- 
lelism offered by the domain decomposition approach 
to solve the coupled fluid-structural interaction prob- 
lems. 


Wing Aeroelasticity 

In order to validate the present development, com- 
putations were done for a clipped delta wing config- 
uration. 27 The transonic flutter characteristics of this 
wing are available from wind tunnel tests for various 
flow parameters. For this computation, the flow field 
is discretized using a C-H grid topology of size 151 x 
30 x 25. The fluid grid is assigned to 32 processors 
on the iPSC/860. The processors are arranged as a 3- 
D mesh of 8 processors in the chordwise direction and 
2 processors in both the spanwise and surface-normal 
directions. 

A 20 degrees-of-freedom (DOF) ANS4 shell 28 ele- 
ment was used for the finite element modeling of the 
structure. The wing is modeled as a plate. Consider- 
ing the wing model used in the experiment, variation of 
mass density is allowed along the chordwise and span- 
wise directions. But the thickness of the finite element 
model is kept constant. This is based on the assump- 
tions that the stiffness of the wing is dominated by an 
aluminum-alloy insert and that the mass distribution of 
the wing is significantly changed due to plastic foams 
covering the aluminum-alloy insert. For the structures 
part of the computation, processors were assigned as a 
2-D mesh of 2 processors in the chordwise and spanwise 
directions, respectively, on the iPSC/860. 

In order to compare sequential and parallel in- 
tegration schemes, the aeroelastic responses were ob- 
tained using both schemes on the iPSC/860. The re- 
sults are presented in Fig. 4. The responses were ob- 
tained for 0° angle of attack (AoA) at = 0.854 and 
a given dynamic pressure of 1.0 psi. The two results 
agree well. For 256 finite elements (1360 DOF.) with 4 
processors on the structural domain, the computational 
times per integration step are 6.22 and 3.33 seconds by 
using sequential and parallel integration schemes, re- 
spectively. A speedup factor of 1.87 is achieved by using 
the parallel integration scheme. The parallel integra- 
tion scheme enables concurrent solution of the coupled 
fluid and structural equations without causing any sig- 
nificant inaccuracy or instability problems. The com- 
putation time per integration step required is deter- 
mined by the computational domain that requires most 
time per integration step. This parallel integration is 
one of the advantages of using MIMD computers for 
multidisciplinary analysis. 

Aeroelastic responses were also computed for var- 
ious other dynamic pressures in order to predict the 
flutter dynamic pressure and compare with the exper- 
iment. Figure 5 shows the stable, near neutrally sta- 
ble, and unstable responses of wing tip displacements 
at the leading edge for dynamic pressures of 0.80, 0.85, 


and 0.90 psi t respectively. From the responses shown in 
Fig. 5, the interpolated dynamic pressure for the neu- 
trally stable condition is 0.84 psi. It is noted that the 
experimental dynamic pressure measured at the neu- 
trally stable condition was 0.91 psi. 27 Considering the 
lack of experimental pressure data on the wing and the 
error involved in modeling the wing as a plate with con- 
stant thickness, the computational result is deemed an 
acceptable prediction of the flutter dynamic pressure. 

Wing-Body Aeroelasticity 

The main purpose of this work is to compute aeroe- 
lastic responses of fully flexible wing-body configura- 
tions on MIMD parallel computers. For this purpose, a 
general-purpose moving-grid capability is required. In 
the present work, an analytical scheme 3 that will gener- 
ate a moving H-0 grid is implemented on the iPSC/860. 
This scheme generates the field grid according to the 
surface grid deformation. For demonstration purposes, 
an HSCT type wing-body configuration (1807 model) 
is selected. Figure 6 show the baseline grid. The size of 
the baseline grid is 95 x 89 x 30. However, it should 
be noted that the technology developed in this work for 
moving grid is independent of grid size. The grid gen- 
erated by the code when the structure is deformed is 
shown in Fig. 7. Note that the singular planes upstream 
of the leading edge and downstream of the trailing edge 
are deformed according to the deformed shape of the 
configuration. 

In order to verify the coupling of the surface move- 
ment with the grid movement, dynamic aeroelastic re- 
sponses are obtained for the above wing-body config- 
uration. Both the body and wing are allowed to be 
flexible. The wing-body configuration is modeled as 
a plate/shell structure using 308 elements. The finite 
element layout is shown in Fig. 8. The structural prop- 
erties are chosen so that the structure can demonstrate 
aeroelastic responses for a given flow condition. Sym- 
metric boundary conditions are applied at the top and 
bottom of the body symmetry lines. All DOF are fixed 
along the bottom symmetry line of the mid-body. This 
results in a total of 1641 DOF for the structure. 

Aeroelastic computations are done on the flexible 
wing-body configuration by directly coupling the pres- 
sures computed solving the Euler equations with the 
FE structural equations of motion. A demonstration 
calculation is done for a dynamic aeroelastic case when 
the configuration is ramping up from 0° to 5° AoA 
at Moo = 2.1 as shown in Fig. 9. This ramping mo- 
tion is started from the steady state of 4.75° AoA and 
Moo — 2.1. It is assumed that the wing root is 300 
inches long and aeroelastic computations are done at 
a dynamic pressure of 3.0 psi. The configuration is 
pitched up about the axis perpendicular to the symme- 


try plane and located at the leading edge of the wing 
root. Starting from the steady state solution, the con- 
figuration is pitched up aba rate of 0.0015 degrees per 
time step. At the end of each time step a new field 
grid is generated that conforms to the deformed sur- 
face. Figure 9 shows the response of the leading edge 
of the tip section. It is noted that the wing continues 
to oscillate after the ramp motion has stopped. This is 
because the inertial force on the structure is still dom- 
inating the aeroelastic motion. 

The effect of aerodynamic forces on the aeroelastic 
responses is studied as shown in Fig. 10. The compu- 
tations are started from the steady state solution at 
4.75° AoA and Moo =2.1 with an initial motion of the 
structure due to uniform accelerations to simulate gust 
loads. When aerodynamic forces are not applied, the 
structure is oscillating without decrease in amplitude. 
But the magnitude of the wing tip deflection is negli- 
gible. Upon applying aerodynamic forces, the wing tip 
deflection becomes significant. However, the aerody- 
namic damping reduces the amplitude of the structural 
responses gradually. 

Performance 

In order to measure the performance of the struc- 
tural domain on the Intel iPSC/860, the FLOP (float- 
ing-point operations) rate on the iPSC/860 is calcu- 
lated by comparing time per integration step on the 
iPSC/860 to the time on the Y-MP using a single 
processor. Operation counts from the Cray Hardware 
Performance Monitor are used. A single processor 
of the iPSC/860 achieve the Y-MP equivalent of 4.2 
MFLOPS, while the corresponding rate is about 77 
MFLOPS on a single Y-MP processor. The Intel rate 
is about 7 percent of the peak performance of a sin- 
gle processor on the iPSC/860. Similar performance 
was reported by Ryan and Weeratunga 21 for the fluid 
domain. All performance data reported are for 64-bit 
arithmetic. 

The performance of the structural domain in par- 
allel ENSAERO has been measured over a wide range 
of processor numbers and problem sizes as shown in 
Fig. 11. The speedup relative to the Y-MP is defined 
as 

t tCray 

speedup = - 

tlntel 

where tcray and t Intel are the computational time 
per integration step measured on the Y-MP and the 
iPSC/860, respectively. Only a single processor is used 
to measure tcray on the Y-MP. The open and filled 
symbols denote the domain decomposition approach 
which results in the minimum and maximum band- 
widths of the stiffness matrix of each subdomain for 
a given number of processors, respectively. 



For the case of 1,360 DOF, the computational time 
per integration step for 64 processors on the iPSC/860 
is close to that on the Y-MP. However, by increasing the 
size of the problem (10,560 and 20,800 DOF), 16 pro- 
cessors of the iPSC/860 achieve about the same speed 
as a single Y-MP processor. It is evident that the JPCG 
solver on the iPSC/860 performs better as the size of 
problem increases. For the case of 20,800 DOF, the rel- 
ative speedup achieved is about 8 when 64 processors 
are in use on the iPSC/860. 

The overall performance of ENSAERO on both the 
Y-MP and the iPSC/860 is shown in Fig. 12 for the case 
of 113,250 grid points for the fluid domain and 10,560 
DOF for the structural domain. In this computation, 
32 processors are assigned to the fluid domain and 16 to 
64 processors to the structural domain. Both the sky- 
line reduction and JPCG solvers are compared on the 
Y-MP while only the JPCG solver is used for the struc- 
tural domain. The height of each column stands for the 
time per integration step. Each col umn is divided into 
times spent for the fluid domain, the structural domain, 
and for idle/intercube communication. 

For the structural domain, it is evident that the 
skyline reduction solver outperforms the JPCG solver 
on the Y-MP. However, the JPCG solver is first imple- 
mented on the iPSC/860. Direct solvers are still un- 
der development on parallel computers. The main pur- 
pose of this work is to compute aeroelastic responses 
of aerospace vehicles on MIMD parallel computers. 
Therefore the well-developed JPCG solver is selected. 
Due to the domain decomposition approach used in this 
work, the JPCG solver can be easily replaced with more 
efficient solvers when they become available. 

When using 32 processors for the structural do- 
main, the JPCG solver on the iPSC/860 achieves the 
performance of the skyline reduction solver on the Y- 
MP. The time per integration step of ENSAERO using 
96 processors in total on the iPSC/860 is about 60% of 
that obtained using the skyline reduction solver with a 
single processor on the Y-MP. This result is based on 
the computation time for the case of 113,250 CFD grid 
points and 10,560 structural equations. It should be 
noted that the structural domain determined the time 
per integration step for this particular problem on the 
iPSC/860. Most of the time on the fluid domain was 
spent waiting for the interface boundary data. How- 
ever, due to the domain decomposition approach, it is 
possible to balance the computational time between the 
two domains by assigning more processors to the struc- 
tural domain. 

Conclusions 

A parallel wing-body version of a multidisciplinary 
code, ENSAERO, has been developed on the Intel 


iPSC/860. A domain decomposition approach was used 
to enable algorithms for the fluid and structural disci- 
plines to be developed and maintained independently. 
This approach provides an efficient and effective envi- 
ronment to researchers. A researcher working in the 
fluid or the structural discipline can develop his own 
algorithms independent of the others. Coupling of the 
disciplines is achieved by exchanging boundary data 
through an intercube communication mechanism. This 
makes it easy for each discipline to incorporate and 
develop new algorithms or data structures without in- 
terferences. 

The performance of the structural domain is far 
behind that of the fluid domain. This is due to the 
less desirable performance of the JPCG algorithm. It 
is noted that direct solvers are still in the early stages of 
development. However, since the procedure developed 
here allows for one domain to select algorithms inde- 
pendent of others, the JPCG algorithm can be easily 
replaced with more efficient algorithms when available. 
Although the solver for the structure is not efficient on 
a serial computer, reasonable computational speed and 
a good load balance can be achieved by assigning more 
processors to the structural domain. The overall time 
per integration step of parallel ENSAERO using 96 pro- 
cessors on the iPSC/860 is reduced to about 60% of the 
best time obtained on a single Y-MP processor for the 
particular problem considered. This shows the advan- 
tage of using the domain decomposition approach for 
the multidisciplinary analysis on MIMD parallel com- 
puters. 

The parallel integration scheme enables the com- 
bination of advanced CFD and CSD technologies with 
minimal increase in computational time per integration 
step while keeping modularity of each discipline. The 
time per integration step is solely determined by the 
domain that requires most computational time on the 
iPSC/860. This parallel integration is one of the advan- 
tages of using MIMD computers for multidisciplinary 
analysis. The procedure developed in this research will 
provide an efficient tool for solving aeroelastic problems 
of complete aerospace vehicle configurations on MIMD 
computers. 
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(a) Fhyoieoi domain (b) Computational domain 


Fig. 1:A schematic diagram of the node-to-element 
approach for the fluid-structural interface. 



(Fluid domain) 
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(a) Sequential Integration 


Fig. 2: Processor arrangements and message ex- 
changes through interprocessor and intercube commu- 
nications. 
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Fig. 4: Aeroelastic responses obtained by using se- 
quential and parallel integration schemes. (M^ = 
0.854, a = 0°, P= 1.0 psi) 



(b) Parallel integration 


Fig. 3: Flow diagrams for sequential and parallel in- 
tegration schemes. 


















Fig. 5: Aeroelastic responses of a clipped delta Fig. 7: Deformed surface and field grids, 

wing obtained by solving finite- difference Euler equa- 
tions and finite-element structural equations of motion. 

(Moo - 0.854, a = 0° ) 



Fig. 6: An HSCT type wing-body configuration with Fig. 8: Finite-element modeling of the wing-body con- 
portions of surface and field physical grids. figuration. 
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Fig. 9: Dynamic aeroelastic displacements at the 
wing-tip leading edge of an HSCT type wing-body with 
ramp motion. (Moo = 2.1, a = 4.75°, P= 3.0 psi) 


Fig. 11: Computational performance of the structural 
domain in ENSAERO with various problem sizes and 
domain decompositions on the Intel iPSC/860 
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Fig. 10: Effect of aerodynamic forces on the aeroelas- Fig. 12: Overall computational performance of EN- 
tic responses of an HSCT type wing-body, (Moo = 2.1, SAERO on the Cray Y-MP and the Intel iPSC/860. 
o = 4.75°) 
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